clc;clear all;close all;
global alpha betta depr sigmaC phi govexp taumax N T ALPHA popweight
global kstst lstst cstst rstst wstst kstart kistart tauK0 mustst
global Delta1 Delta2 lamda transfer kg lam kmax

%mainbenchNoL;
load('benchN2NoL.mat')

%find long-run steady state given initial guesses for lam, Delta1 and Delta2

kistart = [k1ss k2ss];
kstart = sum(popweight.*kistart)+kg;
kmax=5;

%ALPHA=1;

lamda=lam;

Delta2 = 2;
Delta1=-Delta2*phi(2:N)/phi(1);
transfer=0;

mustst = fzero('findststNoL',[-10 10],options)%,0,0)
rstst = 1/betta-1+depr; %r from Euler at the steady state
kstst=e_real/(rstst/alpha)^(1/(1-alpha)); % express k from r=F_k
wstst = (1-alpha)*kstst^alpha*e_real^(-alpha); %w=F_e
lstst=e_real;
cstst = 1/(popweight(1)+sum(popweight(2:N).*lamda))*(kstst^alpha*e_real^(1-alpha) - depr*kstst - govexp); %c from resource constraint
mustst=(cstst^(-sigmaC)+sum(ALPHA.*lamda.*(lamda*cstst).^(-sigmaC)))/(popweight(1)+sum(popweight(2:N).*lamda));

T=200; %we assume that the economy reaches the steady state after T periods

tauK0 = taumax;

%initial guesses 
for i=1:T;
    time1 = i/T;
    X(i) = time1*(kstst) + (1-time1)*kstart; %k
    X(T+i) = mustst; %mu
    X(2*T+i) = max(4-0.1*i,0); %gamma
end;
X(3*T+1:3*T+N) = [Delta1 lamda];
X0=X;

ALPHall=[1.7 1.6 1.65 1.645 1.644 1.6 1.5 1.4 1.3 1.2 1.1 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.05 0 -0.05 -0.1 -0.2 -0.3 -0.32 -0.35 -0.36 -0.37 -0.4 -0.5 -0.7 -0.71 -0.72 -0.73 -0.74 -0.745 -0.76 -0.78 -0.8]
save('solN2NoL.mat')

for ii=1:length(ALPHall)

    ALPHA=ALPHall(ii)
    
    if ii>1
        X=XX(ii-1,:);
    end
    if ii==17 | ii==20 | ii==25
        X=XX(ii-2,:);
    end
    if ii==9
        X=XX(ii-3,:);
    end
    if ii==26
        X=0.5*XX(ii-1,:)+0.5*XX(ii-2,:);
    end
    if ii==33
        X=0.5*XX(ii-1,:)+0.5*X0;
    end
    if ii==34
        X=X0
        for i=1:T
            X(2*T+i) = max(5-0.1*i,0); %gamma
        end
    end
    if ii==35 | ii==36
        X=X0
        for i=1:T
            X(2*T+i) = max(5.2-0.1*i,0); %gamma
        end
    end
    if ii==37
        X=X0
        for i=1:T
            X(2*T+i) = max(5.3-0.1*i,0); %gamma
        end
    end
    if ii==38
        X=X0
        for i=1:T
            X(2*T+i) = max(5.4-0.1*i,0); %gamma
        end
    end
    if ii==39
        X=X0
        for i=1:T
            X(2*T+i) = max(5.5-0.1*i,0); %gamma
        end
    end
    if ii==40
        X=X0
        for i=1:T
            X(2*T+i) = max(6.05-0.1*i,0); %gamma
        end
    end
    if ii==41
        X=X0
        for i=1:T
            X(2*T+i) = max(6.1-0.1*i,0); %gamma
        end
    end
    if ii==42
        X=X0
        for i=1:T
            X(2*T+i) = max(7-0.1*i,0); %gamma
        end
    end 
    if ii==3 | ii==4 | ii==5 | (ii>=10 & ii<=14) | ii==19 | ii==21 | ii==22 | ii==27 | ii==28 | ii==29 | ii==32 | ii==40
        X=0.7*XX(ii-1,:)+0.3*X0;
    end
    if ii==8 | ii==30
        X=0.7*XX(ii-2,:)+0.3*X0;
    end
    if ii==31
        X=0.7*XX(ii-3,:)+0.3*X0;
    end
    Delta1 = X(3*T+1);
    lamda = X(3*T+2:3*T+N);
    if ii<23 | ii>=27% & ii<=34)
        for i=1:T;
            time1 = i/T;
            X(i) = time1*(kstst) + (1-time1)*kstart; %k
        end
    end
    if ii==6
        X=XX(2,:);
    end
Delta2=-Delta1*phi(1)/phi(2:N);
rstst = 1/betta-1+depr; %r from Euler at the steady state
kstst=e_real/(rstst/alpha)^(1/(1-alpha)); % express k from r=F_k
wstst = (1-alpha)*kstst^alpha*e_real^(-alpha); %w=F_e
lstst=e_real;
cstst = 1/(popweight(1)+sum(popweight(2:N).*lamda))*(kstst^alpha*e_real^(1-alpha) - depr*kstst - govexp); %c from resource constraint
mustst=(cstst^(-sigmaC)+sum(ALPHA.*lamda.*(lamda*cstst).^(-sigmaC)))/(popweight(1)+sum(popweight(2:N).*lamda));
    
tauK0 = taumax;
tolf0 = 1e-7;

if ii<=29
    options = optimset('MaxFunEvals',1000000,'MaxIter',100,'TolFun',1e-10,'TolX',1e-10,'Disp','iter');
else
    options = optimset('MaxFunEvals',1000000,'MaxIter',30,'TolFun',1e-10,'TolX',1e-10,'Disp','iter');
end
if ii==34 | ii==35 | ii>=38
    options = optimset('MaxFunEvals',1000000,'MaxIter',40,'TolFun',1e-10,'TolX',1e-10,'Disp','iter');
end

X = fsolve('resid3vNoL',X,options)

    X1=X;
    X=X1;

if ii==9 | (ii>=23 & ii<=26)
    for i=1:T;
        time1 = i/T;
        X(i) = time1*(kstst) + (1-time1)*kstart; %k
    end
    X = fsolve('resid3vNoL',X,options)
end

if ii==2 | ii==12 | ii==17 | ii==18 | ii==19 | ii==21 | ii==22 | ii>=23
    [X,check] = broydn('resid3vNoL',X,tolf0)%,0,1);
    if size(X,1)>1
        X=X';
    end
    sum(resid3vNoL(X).^2)
    X = fsolve('resid3vNoL',X,options)
end

if ii==26 | ii==27 | ii==31 | ii>=34
    [X,check] = broydn('resid3vNoL',X,tolf0)%,0,1);
    if size(X,1)>1
        X=X';
    end
    sum(resid3vNoL(X).^2)
    X = fsolve('resid3vNoL',X,options)
end

if ii==34
     X = fsolve('resid3vNoL',X,options)
end

XX(ii,:)=X;

k = X(1:T);
mu = X(T+1:2*T);
gamma = X(2*T+1:3*T);
Delta1 = X(3*T+1);
lamda = X(3*T+2:3*T+N);
Delta2=-Delta1*phi(1)/phi(2);
rstst = 1/betta-1+depr; %r from Euler at the steady state
kstst=e_real/(rstst/alpha)^(1/(1-alpha)); % express k from r=F_k
wstst = (1-alpha)*kstst^alpha*e_real^(-alpha); %w=F_e
lstst=e_real;
cstst = 1/(popweight(1)+sum(popweight(2:N).*lamda))*(kstst^alpha*e_real^(1-alpha) - depr*kstst - govexp); %c from resource constraint
mustst=(cstst^(-sigmaC)+sum(ALPHA.*lamda.*(lamda*cstst).^(-sigmaC)))/(popweight(1)+sum(popweight(2:N).*lamda));

Delta1*kistart(1)+sum(Delta2.*kistart(2:N))

klast = [kstart X(1:T-1)]; %k_{t-1}
e=ones(1,T)*e_real;
l=e;
lstst=e_real

r = alpha*klast.^(alpha-1).*e.^(1-alpha); %r=F_k
w = (1-alpha)*klast.^(alpha).*e.^(-alpha); %w=F_e
c = 1/(popweight(1)+sum(popweight(2:N).*lamda))*(klast.^alpha.*e.^(1-alpha)+ (1-depr)*klast - k - govexp); %resource constraint

cnext = [c(2:end) cstst];
rnext = [r(2:end) rstst];

%check bound on tauk
constrtest = c.^(-sigmaC)./cnext.^(-sigmaC) - betta*(1+ (rnext-depr)*(1-taumax));

cistst = [cstst lamda*cstst]; %consumption of groups in the long run
lists=[lstst lstst]
for tt=1:T
ci(tt,:) = [c(tt) lamda*c(tt)]; %consumption of groups during the transition
li(tt,:) = [l(tt) l(tt)];
end

%lifetime utilities
if (sigmaC==1)
    V = sum(betta.^[0:T-1]'*ones(1,N).*(log(ci)))+ betta^(T)/(1-betta)*log(cistst);
else
    V = sum(betta.^[0:T-1]'*ones(1,N).*((ci.^(1-sigmaC)-1)/(1-sigmaC)))+ betta^(T)/(1-betta)*(cistst.^(1-sigmaC)-1)/(1-sigmaC);
end

if sum(V>=VsqNoL)==2
    PI(ii)=1
else
    PI(ii)=0
end
Vall(ii,:)=V
VsqNoL
%tax rates
tauh = 1-(sum(betta.^[0:T-1].*(c.^(-sigmaC).*c)) + betta^T/(1-betta)*cstst^(-sigmaC)*cstst-c(1)^(-sigmaC)*kistart(1)*(1+(r(1)-depr)*(1-tauK0)))/(sum(betta.^[0:T-1].*c.^(-sigmaC)*phi(1).*w.*l)+betta^T/(1-betta)*cstst^(-sigmaC)*phi(1)*wstst*e_real);
tauKt = 1 - ((c(1:end-1)./c(2:end)).^(-sigmaC)/betta-1)./(r(2:end)-depr); %Euler
tauKt = [tauK0 tauKt];
tauLt = ones(1,T)*tauh;
taxes = [tauKt tauLt];
[swi dur]=min((tauKt-taumax/2).^2);
dur

figure
plot(tauKt)
save('solN2NoL.mat')

end

%solve for some more values of ALPHA

ALPHall(2)=1.675
ii=2
ALPHA=ALPHall(ii)   
X=XX(ii-1,:);
for i=1:T;
    time1 = i/T;
    X(i) = time1*(kstst) + (1-time1)*kstart; %k
end
    
Delta2=-Delta1*phi(1)/phi(2:N);
rstst = 1/betta-1+depr; %r from Euler at the steady state
kstst=e_real/(rstst/alpha)^(1/(1-alpha)); % express k from r=F_k
wstst = (1-alpha)*kstst^alpha*e_real^(-alpha); %w=F_e
lstst=e_real;
cstst = 1/(popweight(1)+sum(popweight(2:N).*lamda))*(kstst^alpha*e_real^(1-alpha) - depr*kstst - govexp); %c from resource constraint
mustst=(cstst^(-sigmaC)+sum(ALPHA.*lamda.*(lamda*cstst).^(-sigmaC)))/(popweight(1)+sum(popweight(2:N).*lamda));
    
tauK0 = taumax;
tolf0 = 1e-7;

options = optimset('MaxFunEvals',1000000,'MaxIter',100,'TolFun',1e-10,'TolX',1e-10,'Disp','iter');
X = fsolve('resid3vNoL',X,options)
  
X1=X;
X=X1;

[X,check] = broydn('resid3vNoL',X,tolf0)%,0,1);
if size(X,1)>1
    X=X';
end
sum(resid3vNoL(X).^2)
X = fsolve('resid3vNoL',X,options)

XX(ii,:)=X;

k = X(1:T);
mu = X(T+1:2*T);
gamma = X(2*T+1:3*T);
Delta1 = X(3*T+1);
lamda = X(3*T+2:3*T+N);
Delta2=-Delta1*phi(1)/phi(2);
rstst = 1/betta-1+depr; %r from Euler at the steady state
kstst=e_real/(rstst/alpha)^(1/(1-alpha)); % express k from r=F_k
wstst = (1-alpha)*kstst^alpha*e_real^(-alpha); %w=F_e
lstst=e_real;
cstst = 1/(popweight(1)+sum(popweight(2:N).*lamda))*(kstst^alpha*e_real^(1-alpha) - depr*kstst - govexp); %c from resource constraint
mustst=(cstst^(-sigmaC)+sum(ALPHA.*lamda.*(lamda*cstst).^(-sigmaC)))/(popweight(1)+sum(popweight(2:N).*lamda));

Delta1*kistart(1)+sum(Delta2.*kistart(2:N))

klast = [kstart X(1:T-1)]; %k_{t-1}
e=ones(1,T)*e_real;
l=e;
lstst=e_real

r = alpha*klast.^(alpha-1).*e.^(1-alpha); %r=F_k
w = (1-alpha)*klast.^(alpha).*e.^(-alpha); %w=F_e
c = 1/(popweight(1)+sum(popweight(2:N).*lamda))*(klast.^alpha.*e.^(1-alpha)+ (1-depr)*klast - k - govexp); %resource constraint

cnext = [c(2:end) cstst];
rnext = [r(2:end) rstst];

%check bound on tauk
constrtest = c.^(-sigmaC)./cnext.^(-sigmaC) - betta*(1+ (rnext-depr)*(1-taumax));

cistst = [cstst lamda*cstst]; %consumption of groups in the long run
lists=[lstst lstst]
for tt=1:T
ci(tt,:) = [c(tt) lamda*c(tt)]; %consumption of groups during the transition
li(tt,:) = [l(tt) l(tt)];
end

%lifetime utilities
if (sigmaC==1)
    V = sum(betta.^[0:T-1]'*ones(1,N).*(log(ci)))+ betta^(T)/(1-betta)*log(cistst);
else
    V = sum(betta.^[0:T-1]'*ones(1,N).*((ci.^(1-sigmaC)-1)/(1-sigmaC)))+ betta^(T)/(1-betta)*(cistst.^(1-sigmaC)-1)/(1-sigmaC);
end

if sum(V>=VsqNoL)==2
    PI(ii)=1
else
    PI(ii)=0
end
Vall(ii,:)=V
VsqNoL
%tax rates
tauh = 1-(sum(betta.^[0:T-1].*(c.^(-sigmaC).*c)) + betta^T/(1-betta)*cstst^(-sigmaC)*cstst-c(1)^(-sigmaC)*kistart(1)*(1+(r(1)-depr)*(1-tauK0)))/(sum(betta.^[0:T-1].*c.^(-sigmaC)*phi(1).*w.*l)+betta^T/(1-betta)*cstst^(-sigmaC)*phi(1)*wstst*e_real);
tauKt = 1 - ((c(1:end-1)./c(2:end)).^(-sigmaC)/betta-1)./(r(2:end)-depr); %Euler
tauKt = [tauK0 tauKt];
tauLt = ones(1,T)*tauh;
taxes = [tauKt tauLt];

figure
plot(tauKt)

ALPHA=1.8
    
%X=XX(3,:);
X=0.8*XX(1,:)+0.2*XX(2,:);
  
for i=1:T;
    time1 = i/T;
    X(i) = time1*(kstst) + (1-time1)*kstart; %k
end
    
Delta2=-Delta1*phi(1)/phi(2:N);
rstst = 1/betta-1+depr; %r from Euler at the steady state
kstst=e_real/(rstst/alpha)^(1/(1-alpha)); % express k from r=F_k
wstst = (1-alpha)*kstst^alpha*e_real^(-alpha); %w=F_e
lstst=e_real;
cstst = 1/(popweight(1)+sum(popweight(2:N).*lamda))*(kstst^alpha*e_real^(1-alpha) - depr*kstst - govexp); %c from resource constraint
mustst=(cstst^(-sigmaC)+sum(ALPHA.*lamda.*(lamda*cstst).^(-sigmaC)))/(popweight(1)+sum(popweight(2:N).*lamda));
    
tauK0 = taumax;
tolf0 = 1e-7;

options = optimset('MaxFunEvals',1000000,'MaxIter',100,'TolFun',1e-10,'TolX',1e-10,'Disp','iter');
X = fsolve('resid3vNoL',X,options)
  
k = X(1:T);
mu = X(T+1:2*T);
gamma = X(2*T+1:3*T);
Delta1 = X(3*T+1);
lamda = X(3*T+2:3*T+N);
Delta2=-Delta1*phi(1)/phi(2);
rstst = 1/betta-1+depr; %r from Euler at the steady state
kstst=e_real/(rstst/alpha)^(1/(1-alpha)); % express k from r=F_k
wstst = (1-alpha)*kstst^alpha*e_real^(-alpha); %w=F_e
lstst=e_real;
cstst = 1/(popweight(1)+sum(popweight(2:N).*lamda))*(kstst^alpha*e_real^(1-alpha) - depr*kstst - govexp); %c from resource constraint
mustst=(cstst^(-sigmaC)+sum(ALPHA.*lamda.*(lamda*cstst).^(-sigmaC)))/(popweight(1)+sum(popweight(2:N).*lamda));

Delta1*kistart(1)+sum(Delta2.*kistart(2:N))

klast = [kstart X(1:T-1)]; %k_{t-1}
e=ones(1,T)*e_real;
l=e;
lstst=e_real

r = alpha*klast.^(alpha-1).*e.^(1-alpha); %r=F_k
w = (1-alpha)*klast.^(alpha).*e.^(-alpha); %w=F_e
c = 1/(popweight(1)+sum(popweight(2:N).*lamda))*(klast.^alpha.*e.^(1-alpha)+ (1-depr)*klast - k - govexp); %resource constraint

cnext = [c(2:end) cstst];
rnext = [r(2:end) rstst];

%check bound on tauk
constrtest = c.^(-sigmaC)./cnext.^(-sigmaC) - betta*(1+ (rnext-depr)*(1-taumax));

cistst = [cstst lamda*cstst]; %consumption of groups in the long run
lists=[lstst lstst]
for tt=1:T
ci(tt,:) = [c(tt) lamda*c(tt)]; %consumption of groups during the transition
li(tt,:) = [l(tt) l(tt)];
end

%lifetime utilities
if (sigmaC==1)
    V = sum(betta.^[0:T-1]'*ones(1,N).*(log(ci)))+ betta^(T)/(1-betta)*log(cistst);
else
    V = sum(betta.^[0:T-1]'*ones(1,N).*((ci.^(1-sigmaC)-1)/(1-sigmaC)))+ betta^(T)/(1-betta)*(cistst.^(1-sigmaC)-1)/(1-sigmaC);
end

%tax rates
tauh = 1-(sum(betta.^[0:T-1].*(c.^(-sigmaC).*c)) + betta^T/(1-betta)*cstst^(-sigmaC)*cstst-c(1)^(-sigmaC)*kistart(1)*(1+(r(1)-depr)*(1-tauK0)))/(sum(betta.^[0:T-1].*c.^(-sigmaC)*phi(1).*w.*l)+betta^T/(1-betta)*cstst^(-sigmaC)*phi(1)*wstst*e_real);
tauKt = 1 - ((c(1:end-1)./c(2:end)).^(-sigmaC)/betta-1)./(r(2:end)-depr); %Euler
tauKt = [tauK0 tauKt];
tauLt = ones(1,T)*tauh;
taxes = [tauKt tauLt];

figure
plot(tauKt)

term1=(1-sigmaC)*(1-betta)*V(2);
epscap0a=((term1+1).^(1/(1-sigmaC))/cistst_bm(2)-1)*100
term2=(1-sigmaC)*(1-betta)*V(1);
epswor0a=((term2+1).^(1/(1-sigmaC))/cistst_bm(1)-1)*100

ALPHA=1.9

X=0.8*X+0.2*XX(1,:);

for i=1:T;
    time1 = i/T;
    X(i) = time1*(kstst) + (1-time1)*kstart; %k
end
    
Delta2=-Delta1*phi(1)/phi(2:N);
rstst = 1/betta-1+depr; %r from Euler at the steady state
kstst=e_real/(rstst/alpha)^(1/(1-alpha)); % express k from r=F_k
wstst = (1-alpha)*kstst^alpha*e_real^(-alpha); %w=F_e
lstst=e_real;
cstst = 1/(popweight(1)+sum(popweight(2:N).*lamda))*(kstst^alpha*e_real^(1-alpha) - depr*kstst - govexp); %c from resource constraint
mustst=(cstst^(-sigmaC)+sum(ALPHA.*lamda.*(lamda*cstst).^(-sigmaC)))/(popweight(1)+sum(popweight(2:N).*lamda));
    
tauK0 = taumax;
tolf0 = 1e-7;

options = optimset('MaxFunEvals',1000000,'MaxIter',100,'TolFun',1e-10,'TolX',1e-10,'Disp','iter');
X = fsolve('resid3vNoL',X,options)

k = X(1:T);
mu = X(T+1:2*T);
gamma = X(2*T+1:3*T);
Delta1 = X(3*T+1);
lamda = X(3*T+2:3*T+N);
Delta2=-Delta1*phi(1)/phi(2);
rstst = 1/betta-1+depr; %r from Euler at the steady state
kstst=e_real/(rstst/alpha)^(1/(1-alpha)); % express k from r=F_k
wstst = (1-alpha)*kstst^alpha*e_real^(-alpha); %w=F_e
lstst=e_real;
cstst = 1/(popweight(1)+sum(popweight(2:N).*lamda))*(kstst^alpha*e_real^(1-alpha) - depr*kstst - govexp); %c from resource constraint
mustst=(cstst^(-sigmaC)+sum(ALPHA.*lamda.*(lamda*cstst).^(-sigmaC)))/(popweight(1)+sum(popweight(2:N).*lamda));

Delta1*kistart(1)+sum(Delta2.*kistart(2:N))

klast = [kstart X(1:T-1)]; %k_{t-1}
e=ones(1,T)*e_real;
l=e;
lstst=e_real

r = alpha*klast.^(alpha-1).*e.^(1-alpha); %r=F_k
w = (1-alpha)*klast.^(alpha).*e.^(-alpha); %w=F_e
c = 1/(popweight(1)+sum(popweight(2:N).*lamda))*(klast.^alpha.*e.^(1-alpha)+ (1-depr)*klast - k - govexp); %resource constraint

cnext = [c(2:end) cstst];
rnext = [r(2:end) rstst];

%check bound on tauk
constrtest = c.^(-sigmaC)./cnext.^(-sigmaC) - betta*(1+ (rnext-depr)*(1-taumax));

cistst = [cstst lamda*cstst]; %consumption of groups in the long run
lists=[lstst lstst]
for tt=1:T
ci(tt,:) = [c(tt) lamda*c(tt)]; %consumption of groups during the transition
li(tt,:) = [l(tt) l(tt)];
end

%lifetime utilities
if (sigmaC==1)
    V = sum(betta.^[0:T-1]'*ones(1,N).*(log(ci)))+ betta^(T)/(1-betta)*log(cistst);
else
    V = sum(betta.^[0:T-1]'*ones(1,N).*((ci.^(1-sigmaC)-1)/(1-sigmaC)))+ betta^(T)/(1-betta)*(cistst.^(1-sigmaC)-1)/(1-sigmaC);
end

%tax rates
tauh = 1-(sum(betta.^[0:T-1].*(c.^(-sigmaC).*c)) + betta^T/(1-betta)*cstst^(-sigmaC)*cstst-c(1)^(-sigmaC)*kistart(1)*(1+(r(1)-depr)*(1-tauK0)))/(sum(betta.^[0:T-1].*c.^(-sigmaC)*phi(1).*w.*l)+betta^T/(1-betta)*cstst^(-sigmaC)*phi(1)*wstst*e_real);
tauKt = 1 - ((c(1:end-1)./c(2:end)).^(-sigmaC)/betta-1)./(r(2:end)-depr); %Euler
tauKt = [tauK0 tauKt];
tauLt = ones(1,T)*tauh;
taxes = [tauKt tauLt];

figure
plot(tauKt)

term1=(1-sigmaC)*(1-betta)*V(2);
epscap0b=((term1+1).^(1/(1-sigmaC))/cistst_bm(2)-1)*100
term2=(1-sigmaC)*(1-betta)*V(1);
epswor0b=((term2+1).^(1/(1-sigmaC))/cistst_bm(1)-1)*100

ALPHA=2

X=0.8*X+0.2*XX(1,:);

for i=1:T;
    time1 = i/T;
    X(i) = time1*(kstst) + (1-time1)*kstart; %k
end
    
Delta2=-Delta1*phi(1)/phi(2:N);
rstst = 1/betta-1+depr; %r from Euler at the steady state
kstst=e_real/(rstst/alpha)^(1/(1-alpha)); % express k from r=F_k
wstst = (1-alpha)*kstst^alpha*e_real^(-alpha); %w=F_e
lstst=e_real;
cstst = 1/(popweight(1)+sum(popweight(2:N).*lamda))*(kstst^alpha*e_real^(1-alpha) - depr*kstst - govexp); %c from resource constraint
mustst=(cstst^(-sigmaC)+sum(ALPHA.*lamda.*(lamda*cstst).^(-sigmaC)))/(popweight(1)+sum(popweight(2:N).*lamda));
    
tauK0 = taumax;
tolf0 = 1e-7;

options = optimset('MaxFunEvals',1000000,'MaxIter',100,'TolFun',1e-10,'TolX',1e-10,'Disp','iter');
X = fsolve('resid3vNoL',X,options)

[X,check] = broydn('resid3vNoL',X,tolf0)%,0,1);
if size(X,1)>1
    X=X';
end
sum(resid3vNoL(X).^2)

k = X(1:T);
mu = X(T+1:2*T);
gamma = X(2*T+1:3*T);
Delta1 = X(3*T+1);
lamda = X(3*T+2:3*T+N);
Delta2=-Delta1*phi(1)/phi(2);
rstst = 1/betta-1+depr; %r from Euler at the steady state
kstst=e_real/(rstst/alpha)^(1/(1-alpha)); % express k from r=F_k
wstst = (1-alpha)*kstst^alpha*e_real^(-alpha); %w=F_e
lstst=e_real;
cstst = 1/(popweight(1)+sum(popweight(2:N).*lamda))*(kstst^alpha*e_real^(1-alpha) - depr*kstst - govexp); %c from resource constraint
mustst=(cstst^(-sigmaC)+sum(ALPHA.*lamda.*(lamda*cstst).^(-sigmaC)))/(popweight(1)+sum(popweight(2:N).*lamda));

Delta1*kistart(1)+sum(Delta2.*kistart(2:N))

klast = [kstart X(1:T-1)]; %k_{t-1}
e=ones(1,T)*e_real;
l=e;
lstst=e_real

r = alpha*klast.^(alpha-1).*e.^(1-alpha); %r=F_k
w = (1-alpha)*klast.^(alpha).*e.^(-alpha); %w=F_e
c = 1/(popweight(1)+sum(popweight(2:N).*lamda))*(klast.^alpha.*e.^(1-alpha)+ (1-depr)*klast - k - govexp); %resource constraint

cnext = [c(2:end) cstst];
rnext = [r(2:end) rstst];

%check bound on tauk
constrtest = c.^(-sigmaC)./cnext.^(-sigmaC) - betta*(1+ (rnext-depr)*(1-taumax));

cistst = [cstst lamda*cstst]; %consumption of groups in the long run
lists=[lstst lstst]
for tt=1:T
ci(tt,:) = [c(tt) lamda*c(tt)]; %consumption of groups during the transition
li(tt,:) = [l(tt) l(tt)];
end

%lifetime utilities
if (sigmaC==1)
    V = sum(betta.^[0:T-1]'*ones(1,N).*(log(ci)))+ betta^(T)/(1-betta)*log(cistst);
else
    V = sum(betta.^[0:T-1]'*ones(1,N).*((ci.^(1-sigmaC)-1)/(1-sigmaC)))+ betta^(T)/(1-betta)*(cistst.^(1-sigmaC)-1)/(1-sigmaC);
end

%tax rates
tauh = 1-(sum(betta.^[0:T-1].*(c.^(-sigmaC).*c)) + betta^T/(1-betta)*cstst^(-sigmaC)*cstst-c(1)^(-sigmaC)*kistart(1)*(1+(r(1)-depr)*(1-tauK0)))/(sum(betta.^[0:T-1].*c.^(-sigmaC)*phi(1).*w.*l)+betta^T/(1-betta)*cstst^(-sigmaC)*phi(1)*wstst*e_real);
tauKt = 1 - ((c(1:end-1)./c(2:end)).^(-sigmaC)/betta-1)./(r(2:end)-depr); %Euler
tauKt = [tauK0 tauKt];
tauLt = ones(1,T)*tauh;
taxes = [tauKt tauLt];

figure
plot(tauKt)

term1=(1-sigmaC)*(1-betta)*V(2);
epscap0c=((term1+1).^(1/(1-sigmaC))/cistst_bm(2)-1)*100
term2=(1-sigmaC)*(1-betta)*V(1);
epswor0c=((term2+1).^(1/(1-sigmaC))/cistst_bm(1)-1)*100

save('solN2NoL.mat')

clc;clear all;close all;
load('solN2NoL_fb.mat')
load('solN2NoL.mat')

PI
Vall
VsqNoL
Vcap=Vall(:,2)
Vwor=Vall(:,1)
 
term1=(1-sigmaC)*(1-betta)*Vcap;
epscap=((term1+1).^(1/(1-sigmaC))/cistst_bm(2)-1)*100
term2=(1-sigmaC)*(1-betta)*Vwor;
epswor=((term2+1).^(1/(1-sigmaC))/cistst_bm(1)-1)*100
epscap=[epscap0c epscap0b epscap0a epscap']
epswor=[epswor0c epswor0b epswor0a epswor']


epscap
epswor

plot(epswor(8:26),epscap(8:26),'-k','LineWidth',2.5);
axis([-0.2 6 -0.2 3.3])
hold on
plot(epswor(1:7),epscap(1:7),'-.k','LineWidth',2);
hold on
plot(epswor_fb(4:18),epscap_fb(4:18),'--k','LineWidth',1.5);
hold on
plot(epswor_fb(1:4),epscap_fb(1:4),'-.k','LineWidth',1);
hold on
plot(epswor(14),epscap(14),'-k','Marker','.','MarkerSize',25,'LineWidth',1.5);
hold on
plot(epswor(26),epscap(26),'-k','Marker','s','MarkerSize',7,'MarkerFaceColor','k','LineWidth',1.5);
hold on
plot(epswor_fb(19:21),epscap_fb(19:21),'-.k','LineWidth',1);
hold on
plot(epswor(27:34),epscap(27:34),'-.k','LineWidth',2);
hold on
plot([-0.2 6],[0 0],'--k');
hold on
plot([0 0],[-0.2 3.3],'--k');
xlabel("workers' welfare increase (percent)",'FontSize',12)
ylabel("capitalists' welfare increase (percent)",'FontSize',12)
legend('POPI','PO, not PI','first-best PI','first-best, not PI')

